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Abstract 

It ha been proved that the lattice Boltzmann Methods (LBMs) are very efficient for computational fluid dynamics 
purposes, and for capturing the dynamics of weak acoustic fluctuations. Especially, the multi-relaxation-time lattice 
Boltzmann method (MRT-LBM) appears as a very robust scheme with high precision. It is well-known that there 
exist several free relaxation parameters in the MRT-LBM. Although these parameters have been tuned via linear 
analysis, the sensitivity analysis of these parameters and other related parameters are still not sufficient for detecting 
the behaviors of the dispersion and dissipation relations of the MRT-LBM. Previous researches have shown that the 
bulk dissipation in the MRT-LBM induces a significant over-damping of acoustic disturbances. This indicates that 
MRT-LBM cannot be used to obtain the correct behavior of pressure fluctuations because of the fixed bulk relaxation 
parameter. In order to cure this problem, an effective algorithm has been proposed for recovering the linearized 
Navier-Stokes equations (L-NSE ) from the linearized MRT-LBM (L-MRT-LBM). The recovered L-NSE appear as 
in matrix form with arbitrary order of the truncation errors with respect to St. Then, in wave-number space, the 
first/second-order sensitivity analyses of matrix eigenvalues are used to address the sensitivity of the wavenumber 
magnitudes to the dispersion-dissipation relations. By the first-order sensitivity analysis, the numerical behaviors 
of the group velocity of the MRT-LBM are first obtained. Afterwards, the distribution sensitivities of the matrix 
eigenvalues corresponding to the linearized form of the MRT-LBM are investigated in the complex plane. Based on 
the sensitivity analysis and the recovered L-NSE, we propose some simplified optimization strategies to determine 
the free relaxation parameters in the MRT-LBM. Meanwhile, the dispersion and dissipation relations of the optimal 
MRT-LBM are quantitatively compared with the exact dispersion and dissipation relations. At last, some numerical 
validations on classical acoustic benchmark problems are shown to assess the new optimal MRT-LBM. 

Keywords: Computational aeroacoustics, BGK-LBM, MRT-LBM, Sensitivity, Group velocity, Eigenvalue 
distribution, Dispersion, Dissipation, Von Neumann analysis 



1. Introduction 

The lattice Boltzmann method (LBM) has emerged as an innovative mescoscopic numerical method for the com- 
putational modelling of a wide variety of complex fluid flows [1]. The early researches about the dispersion and 
dissipation relations of the classical LBM (BGK-LBM) show that the LBM are well suited to capture the weak acous- 
tic pressure fluctuations [2, 3, 4, 5]. But, the BGK-LBM always suffers from an unstability drawback. In order to 
improve the stability of the BGK-LBM, the MRT-LBM was derived in moment space based on the generalized lattice 
Boltzmann equation [6, 7]. Thanks to its adjustable bulk viscosity, the MRT-LBM can guarantee a better stability by a 
high value of the bulk viscosity [3,7]. In most researches and simulations, the relaxation parameters of the MRT-LBM 
were chosen according to the results of the linear analysis [6, 7] . Results of the linear von Neumann analysis of the 
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BGK-LBM and the MRT-LBM yielded quantitative results about the dispersion and dissipation relations. Von Neu- 
mann analysis shows that the LBM have lower numerical dissipation than many aeroacoustic-optimized high-order 
schemes for the Navier-Stokes equations [3]. The 2nd-order accurate LBM has the better dispersion capabilities than 
the classical Navier-Stokes schemes with the 2nd-order accuracy in space combined with 3-step Runge-Kutta time 
integration [3]. The dispersion errors in the LBM are larger than those induced by 3rd-order accurate in space the 
finite difference method combined with 4th-order accurate time-integration methods, and also higher than those of 
6th-order accurate dispersion relation preserving (DRP) schemes. But it is important to emphasize that thanks their 
compactness, for a given dispersion error, the LBM are faster than high-order Navier-Stokes equations [3] and are 
therefore more efficient for practical computations. Compared with the BGK-LBM, it is observed that there exists a 
higher dissipation of the acoustic modes in the MRT-LBM [3, 8]. The free-parameter sensitivity of the MRT-LBM is 
missing in open literature. Therefore, it is difficult to guarantee that the the results of the linear analysis are always 
acceptable. In order to clarify the sensitivity issue, we will use a novel methodology to reflect the relation between 
the L-NSE and the L-MRT-LBM. This methodology offers us a way to determine the free relaxation parameters in the 
MRT-LBM. This entirely new work is different from the linear analysis. Using the linear analysis, it is impossible to 
obtain the deterministic values of free relaxation parameters [6, 7]. This difficulty is attributed to a drawback of the lin- 
ear analysis of L-MRT-LBM: the linear analysis [7] cannot reflect the influence of the linearized high-order truncated 
terms corresponding to the L-NSE on the dispersion and dissipation relations. But, using our new methodology, the in- 
fluence of linearized high-order terms coupled with free relaxation parameters on dispersion and dissipation relations 
is revealed under a certain accuracy with respect to the higher-order truncation errors. Meanwhile, this methodology 
offers us a way to deal with the acoustic problems for any values of the shear viscosity and the bulk viscosity. The 
important fact is that the optimal MRT-LBM can be used in the weakly-compressible acoustic problems, for which 
the shear viscosity is not always equal to the bulk viscosity. It is well-known that the classical BGK-LBM can only 
be used to solve acoustic problems with equal shear viscosity and bulk viscosity, while the original MRT-LBM uses 
fixed bulk relaxation parameters. As mentioned above, if the compressible effects are considered, the bulk viscosity 
becomes important for the sound wave propagation [9, 10, 11]. So, it is necessary for acoustic problems to specify 
both the shear viscosity and the bulk viscosity in the LBM. 

By our novel methodology, the recovered L-NSE and the exact L-NSE can be expressed in matrix form in the 
wave-number space as follows 

d t W = B {n) • W + 0(5f) (recovered L - NSE), d t W = B • W (exact L - NSE), (1) 

where the square matrices and B are functions of wave-number k and base-flow uniform velocity u, and B^ 
can be regarded as a perturbation of B [8]. The vector W denotes the perturbed macroscopic fluid flow conservative 
quantities (such as, density p, momentum j, • • •)• Denoting = B^ - B the perturbation matrix, the matrix 

will be the function of all relaxation parameters of the MRT-LBM for n > 2. With the aid of the matrix 
B^ n \ we only consider the influence of the free parameters in the L-MRT-LBM on the hydrodynamic modes. This 
is of great benefit to determine the free relaxation parameters by our proposed optimization strategies. In order to 
establish the optimization strategies, the parameter sensitivities are investigated in the wave-number space. Then, 
the sensitivity properties of the free parameters on the hydrodynamic modes and the non-hydrodynamic modes are 
obtained. According to the sensitivity analysis, we propose a class of the simplified optimization strategies to minimize 
the dispersion and dissipation errors for the MRT-LBM in the wave-number space, based on the matrix perturbation 
theory and the theory of the modified equations. The proposed methodology leads to a relation between the recovered 
L-NSE and the exact L-NSE, and free parameters can be determined by this process. Then, the optimal dispersion 
and dissipation relations are investigated by von Neumann analysis. It is noted that this famous and reliable analysis 
method has been revisited and extended [12]. 

It is necessary to point out that our proposed original methodology is more general than the usual method for 
recovering the higher-order lattice Boltzmann schemes [13, 14], although we also use the Taylor expansion method. 
Our method appears as an easy-to-hand recursive algorithm in the wave-number space, which is completely different 
from the method proposed in [13]. 

In next section, the fundamentals of the LBM are surveyed, and the algorithm used to establish the transformation 
relation from the L-MRT-LBM to the L-NSE is given. In the third section, with the aid of von Neumann method, 
the sensitivity properties of the free parameters to the dispersion and dissipation relations are analyzed in the wave- 
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number space. Then, the optimization strategies are studied in section 4. In the last section, the optimal MRT-LBM 
schemes are validated considering some benchmark problems. 



2. Fundamentals of the LBM and the methodology from the L-MRT-LBM to the L-NSE 

In this section, the fundamental theory of the LBM is briefly recalled. Then, we give the linearized LBM and show 
the method to establish the relation between the linearized LBM and the L- NSE. 

2.1. Lattice Bolt zmann schemes 

The governing equations of the lattice Boltzmann schemes are described as follows [6, 7, 8] 

fiix + Vi6t, t + 6t) = Mx, t) + Ay (/f q) (x, t) - fj(x, 0) , < i, j < N, (2) 

where belongs to the discrete velocity set r V, fi(x, t) is the discrete single particle distribution function corre- 
sponding to Vj and ff q) denotes the discrete single particle equilibrium distribution function. St denotes the time step 
and N + 1 is the number of discrete velocities. Ay is the generalized relation matrix. From here on, the repeated index 
indicates the Einstein summation is used except for some special explanations. Let £ e R d (d denotes the spatial 
dimension) denotes the lattice system, and the following condition is required [13] 

x + vjSt e £, (3) 

that is to say, if x is a node of the lattice, x + VjSt is necessarily another node of the lattice. 
Generally, for BGK-LBM, the relaxation matrix is given by 

Ay = s6 ij9 (4) 

where s denotes the relaxation frequency of BGK-LBM. 

The standard isothermal MRT-LBM has the following form [7, 13] 

m i = W i = mf q \0<i<d, (5) 

and 

ntiix + Stvj, t + 6t) = nti(x, t) + s t (mf q \x, t) - m f (jc, i)\,d+\<i<N. (6) 

It is necessary to point out that for isothermal flow, the number of the conservative quantities is equal to d + 1. 
According to the work of Lallemand and Luo [7], the relaxation parameters in Eq. (6), should satisfy the following 
stability constraints 

j f G(0,2),d+l <i<N. (7) 
The relaxation matrix A corresponding to Eqs. (5) and (6) is defined by 

A = M~ l SM, (8) 

where S is a diagonal matrix which denotes the relaxation parameters of MRT-LBM. M = [Mi A is the 

fc r V lJ J0<i<N,0<i<N 

transformation matrix, which satisfies the following basic conditions (refer to Appendix A for the detail definitions) 
[7] 

M j - 1, M aj - v",(l < a < d). (9) 
The macroscopic quantities are denned by [7, 13] 

m i = M ij fj, mf q) = M^. (10) 
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2.2. Theoretical dispersion and dissipation relations for the L-MRT-LBM and the L-NSE 

In order to implement von Neumann analysis for the L-MRT-LBM, the expression of the L-MRT-LBM in the 
frequency-wave number space is given in this part. Considering a uniform mean part ff and a fluctuating part f/, the 
equilibrium distribution function can be linearized as [3, 4] 



f? q \{j?+f/h^N) = f™'° 



df 



(eq) 



dfj 



fi=fj 



■jr +<>(/?). 



Then, considering a plane wave solution of the linearized equation 

f'j = A 7 exp[i(k • x - cot)], 



(11) 



(12) 



by Eqs. (11), (12) and (10), we get the following eigenvalue problem for the L-MRT-LBM in the frequency- wave 
number space [3, 4, 7] 



M mrt f 



where the matrix iT =A~ l [l- M~ l SM and N b % k is defined by 



N ^ = 6ij - df * 



•(eq) 



l J 



fj 



(13) 



(14) 



(15) 



For the L-NSE, the analytical acoustic modes 61/ (k) and shear modes co s (k) are given by [15] 

Rel^fk)] = \k\(±c s + |u|cos(k^u)), 
/ m [^(k)] = -|k| 2 i(^y + 77), 

Re[oj s (k)] = |k||u|cos(k • u), 
Im[co s (k)] = -|k| 2 v, 

where v is the shear viscosity and n is the bulk viscosity. 

2.3. New form of the L- MRT-LBM and corresponding the higher-order L- NSE 

The equilibrium function mf q) (d + 1 < i < N) is expressed as a function of the conservative quantities Wj 
(Wi = m (eq) ,0 < i < d) (we use the same notations as Dubois and Lallemand [13]) 



or 



mf q) = d ({Wj}o<j<d) ,d+\<i<N. 
mf q) = Gt ({Wjh^ d ) = GijWj = G ljmj 



(16) 



(17) 



In order to implement the linear stability analysis and recover linearized macroscopic equations, we introduce the 
linearized form of Eq. (16) around the reference state [3, 7, 13]. Using Eq. (17), the linearized description of Eqs. (5) 
and (6) can be written as follows 



ntiix, t + 6t) = lx P ' pr m r (x - viSt, t), 



where the matrix has the following form 



V lJ /0<i<d,0<j<d 

® {iv-s tJ ) d 



(18) 



(19) 



/d+l<i<N,d+l<j<N ) 

where @/ 7 = S[Gij (A matrix with the size (N - d - 1) x (d + 1)), I is an identity matrix and the diagonal matrix S is 
defined by 

S = Diag(a_0, s d+u ...,s N ). (20) 

d+l N-d-l 
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In order to derive the linearized high-order equations, one assumes that the discrete single particle distribution 
belongs to C°°(T x £) (a functional set, in which the element possesses a sufficiently smooth property with respect to 
the time domain T and the spatial domain £ ). This assumption is also used by Junk et. al. for asymptotic analysis of 
LBM [16]. This regularity hypothesis indicates that macroscopic quantities are smooth ones and that the linearized 
system (18) is well defined. 

The next step consists of performing Taylor series expansion of the right hand side of Eq. (18), yielding 



v-i St n 

(x, t + St) = J] —Mui-vydaTM^prmr, (21) 



where a indicates the summation from 1 to d. 

Now, we define the matrix A (w) '* = (A^)'* ) as follows 

V i J /0<i<N,0<j<N 



Af'* = ^M a {-v<!d a ) n MJ^ pr . (22) 

When we need to derive the equivalent equations or the modified equations, it is difficult to use the matrix A (w) '* to 
carry out the calculations. In order to overcome this difficulty, we use the differential operators in spectral space. Let 
us note d a = ik a , with k a the wave-number in the a -direction and i 2 = -1. Then, in spectral space, the matrix A (w) '* 

has the following form (a (h) = (A (n) ) ) 

fc \ V i J /0<i<N,0<j<N) 

^^.''(-iftf^V (23) 
Therefore, Eq. (21) can be rewritten as follows 

j-i 

rrii(x, t + 6t)=Y St n A (n) m r + 0(St J ). (24) 



n=0 



In order to derive L-NSE corresponding to L-MRT-LBM defined by Eq. (24), we introduce an original recursive 
algorithm. Given m ; = W/(0 < i < d) (W; denotes the macroscopic conservative quantities, d indicates the number of 
the macroscopic conservative quantities), the algorithm is given as follows 



Initial step. The initial (1) and are given as follows 



<D^ = Sij(0 < i < d\ = -Wij(d + 1 < i < N), (25) 
B W =A ( W. (26) 

ij ir rj v > 

Let W = {Wi)o<i<d and m = {m/}o<i<iv denote the vector of the conservative quantities and the vector of all 
macroscopic quantities, respectively. 

At the first order of St, for all macroscopic quantities, we have 

mi = Q>\fWj + 0(St), < i < N. (27) 

In the matrix form, we obtain 

m = ® (1) W + 0(St). (28) 
At the first-order of St, for conservative quantities, we have 

d t Wi=A^^Wj + 0{St). (29) 

The matrix form is 

d t W = A (l) ® (l) Wj + 0(Si). (30) 
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• Recursive formula for all macroscopic quantities. can be given as follows 

i j 



o<?> = - 



( n-l n-l \ 

01 _ / 



,</+l<i<N (31) 



<D<"?=5 y ,(0< «<<*). (32) 



/=1 " /: 

and 

Eliminating the higher-order term of St n ~ l , we have 

n-l 
Z=0 

where Coeff(-, ■, •) is a function which extracts the coefficients of the polynomials, for example, f(x) = Y!i=o a i %l 

Coeff(/(jc),jc,0 = Of. (34) 

According to Eqs.(31) and (33), we have 

m = (w) W+6>(£f). (35) 
• Recursive formula for conservative quantities. B*f) is presented as follows 



2#> = - V t-*- + V 5t l - l A^ n+l -\ (36) 

Eliminating the higher-order term of c)^ -1 , we have 

n-l 

Efg = Yj ^Coeff (B\f, St, l) . (37) 

1=0 

Now, for the conservative quantities, we have the following equation system 

d t W = B {n) • W + 0(St n ). (38) 

Using Eqs. (33), (35), (37) and (38), we can get the coefficient matrix of the conservative quantities at any order 
of St . Details and validations are displayed in [8]. 

2.4. Illustrative examples 

In this section, we will use the above proposed algorithm to recover the L-NSE. Here, we will show two examples 
corresponding to 2D and 3D MRT-LBM that are similar to the incompressible NSE [6, 7] . 

2.4.7. Application to the 2D MRT-LBM 

In this part, we illustrate the algorithm presented in Sec. 2.3 considering the 2D MRT-LBM. For the standard 2D 
MRT-LBM, the equilibrium distribution functions are described as follow (similary to the incompressible LBS) [7] 



m 



(eq) 



= { P , j X , jy, -2 P + 3 (f, + f y ) , P - 3 (f x + f y ) , -j X , ~jy, (f x - f y ) , j X jy) , (39) 



where j x and j y denote the x-momentum and y-momentum respectively, and p represents the density (the reference 
density po = h Wq - mo = p, W\ - m\ - j x = u, W2 = m<i - j y = v). The corresponding matrix ¥ is given in 
Appendix A. 1 . The diagonal elements of the corresponding diagonal matrix S are set as follows 



S0 = Si = S 2 = 0, S3 = S e , S 4 = S € , S 5 = S 6 = S q , S 7 = S% = S v . 
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(40) 



For the original MRT-LBM [7], only s v is a free parameter, and s e = 1.64, s £ = 1.54, s q = 1.9. The analogous form of 
*F can be found in the existent literature [7, 13]. The derivation of ¥ can be achieved by means of the first-order Taylor 
series expansion with respect to p, j x and j y at reference states [3, 7, 8]. In the expression of X F, (U, V) to denote the 
mean flow velocity components. In order to express the L-NSE by intuitive parameters, we introduce the following 
relation 

^ = v n -i (41) 

where 77 stands for any notations in the set S = {e, e, q, v}. Now, when the truncated error term is equal to 0(St 2 ) 
, the coefficient matrix with the mean flow can be described by the summation of two matrices (the coefficient 
matrices of S° and St) [8]. These two matrices describe the specific terms in the linearized Navier-Stokes equations. 
The coefficient matrix associated with St° is 







-ik x 



-\ik x -2ik x U-ik y V 



-ik v 



-ik y U 



3 1K y 



-ik x V -2ik v V-ik x U 



(42) 



The coefficient matrix associated with St is given in Appendix B.l. The higher-order truncated error matrices of St 
are omited for the sake of brevity (see [8] for more detailed descriptions). 

2.4.2. Application to the 3D MRT-LBM 

In this part, the algorithm presented in Sec. 2. 3 will be used to recover the 3D L-NSE. For the sake of convenience, 
we consider the standard 3D MRT-LBM and the equilibrium distribution functions with 15 discrete velocities are 
given by [6, 7] 



in 



(eq) 



= jp, jx, jy, jz, ~P + {fx + jy + jfj ' ~P> ~\jx> ~\jy> ~\jz-> 2 jl ~ jy ~ f z * fy ~ f z > jxjy, jyjz, jxjz, ^ j • 



(43) 



where j x , j y and j z denote the x-momentum and y-momentum respectively, and p represents the density (the reference 
density po = 1, Wo = mo = p, W\ = m\ = j x = u, W2 = rri2 = j y = v, W3 = = j z = w). The relaxation parameters 
corresponding to the matrix (20) is defined by 



{0, 0, 0, 0, S4, s 5 , S 6 , S 7 , Ss, Sg, Sw, Sn, S12, S13, Su] 

where 

^4 = Se, $5 = S e , S 6 = S 7 = S% = S q , Sg - S\() = S\\ 

For the original MRT-LBM [6, 7], s v is a free paramter and s e = 1.6, s ( 



(44) 



Sl2 = S13 

Z 1.2, Sn 



s t . (45) 

q - 1.6 , s t = 1.2 [6]. For convenience, 
we redefined an enlarged set S = {e, 6, q, v, t). The cooresponding defination (41) can also be used for 3D MRT-LBM. 
When the truncated error term is equal to 0(St 2 ), the coefficient matrix B^ with the mean flow can be expressed by 
the summation of two matrices. The coefficient matrix corresponding to St° is 





-1/3 ik x 

-1/3 iky 

-1/3 ik z 



-2ik x U-ik 7 W- 



-ik x 
ikyV 
-ik x V 
-ik x W 



-ik v 



-\k 7 



-ikyU 



-2ik v V-ik x U -ik 7 W 



-ik z U 
-ik z V 

-iky W -ik x U - ik y V - 2ik z W 



(46) 



The coefficient matrix corresponding to St is given in Appendix B.2 . 

In Sec. 2.4.1 and Sec 2.4.2, the recovered L-NSE are given and the corresponding truncated error is up to St 2 . In 
order to save space, the matrices B^ with the higher-order truncated terms of St n (n > 2) are not given in this paper. 
For readers, by the algorithm given in Sec. 2.3 and Maple symbolic calculations, it is easy to obtain the matrix B^ 
with any value of n. 

From the matrices corresponding to the dissipation, it is clear that the dissipation action of the recovered L-NSE 
is dependent on the mean flow. This dependency always leads to the breakdown of Galilean invariance. This point 
remains an open issue to the knowledge of the authors. 
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3. Sensitivity and stability analysis of free parameters in MRT-LBM 



In this section, the sensitivity of the free parameters and the stability for the MRT-LBM will be addressed by 
means of the first/second-order sensitivity of the eigenvalues [17] and the distribution of the matrix eigenvalues in the 
complex plane, respectively. The free relaxation parameters in E determined by the linear analysis in the wave-number 
spaces are "sub-optimal" [6, 7]. However, the original sub-optimal free relaxation parameters are still used for the 
simulations of fluid flows in open literature. Because of the linear character of the analysis, the influence of the free 
relaxation parameters on the hydrodynamic modes is still not completely clarified. In this part, all of these will be 
analyzed. The analysis offers us a new methodology to propose our simplified optimization strategies for determining 
the free relaxation parameters. 

3.1. Methodology for sensitivity and stability analysis 

In this part, the methods for investigating the sensitivity and stability are presented. 

3.1.1. The fir st/ second-order derivatives of matrix engenvalues with respect to free parameters 

For the convenience of analysis, let T denotes the free parameter set for both 2D and 3D problems. From Eq. (13), 
the linear behavior of MRT-LBM is fully determined by the eigvalues of the matrix M™ 1 . Let X M = {xi}o<i<N and 
Am = {^i)o<i<N denote the eigenvector set and the eigenvector set of the matrix M 1 ™" 1 , respectively. It is clear that X M 
and Am are function of T. Let s e T denote the the concerned parameter which exists in the matrix M 1111 " 1 . Because 
the eigenvalues of the matrix M™^ are distinct from each other (|k| > 0), the first-order sensitivity method can be used 
to address the linear response behaviors or the group velocity behaviors. According to the definition of the eigenvalue 
problem, the following equations are satisfied 

M mrt (£)x/(£) = MsMs). (47) 

Differentiating with repect to s yields 

M^&mie) + M^(s)x[(s) = %(e)Xi(s) + Us^s), (48) 

where Mf^is) = dgM^is). If we define Y M = X M * = {yf }o<i<j+i ifl denotes the conjugate transpose) as the left 
eigenvector set corresponding to A M , we obtain 

Y^(s)M™Xs)X M (s) + Y^(s)M mn (s)X' M (s) = A' M (e) + Y^(s)X'(s)A M (s). (49) 
Considering the diagonal elements in Eq. (49), we have 

yf (s)M™Xs)xi(s) + yf (s)M mrt (s)x^s) = A[(s)yf (s)x j {s) + Ms)yf ( e )xJ(e). (50) 

Since yfx t = 1, we get 

A\{s) = yf {s)Mf(s)x i {s). (51) 

It is clear that the derivatives of eigenvalues with respect to the free parameters can be expressed by the right and 
left eigenvectors and the derivative of the matrix M 1 ™ 1 . The first-order derivative of eigenvectors can be described by 
[18, 19] 

x\(s) = Qjxjis), yf '(g) = -Cjijjis) (52) 

where the matix C is defined by 

c f yf(s)Mf It (s)x i (s)/(A i (s)-A j (s)), itj 
11 1 0, i - j 

From Eq. (51), the second-order derivative can be described by [18, 19] 

A'iis) = yf (£)M™ rt (e) Xl (£) + yf (e^sWe) + yf(e)M™ rt (£)x;(e). (54) 
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From Eq. (13), the following relation between co = {k>j}o<i<iv and A M holds 

exp(-ioji(s)) = Ai(s). (55) 

So, the first and second order derivatives of a)i(s)) can be expressed by 

0)[{s) = iX^IUsX oj'/(s) = i^(s)/Ai(s) - i^isf/Ms) 2 . (56) 

Now, the first/second-order sensitivity of free parameters to 6^(-) can be investigated easily. The behavior of the 
numerical group velocity of the MRT-LBM can be addressed. 



3.1.2. Hydro dynamic stability in the complex planes 

It is known that the changes of any free parameters always result in the variation of the matrix M™ 1 . By inves- 
tigating the eigenvalue distribution of the matrix M™* in the complex planes, the basic sensitivity behaviors of the 
L-MRT-LBM can be revealed. According to Eq. (13), the frequency a> and the eigenvalues of the matrix M™* 1 satisfy 
the following relation 

a)i = i- Log(Ai). (57) 

The stability of MRT-LBM requires that the imaginary part of a>i should be smaller than or equal to 0. It is well-known 
that Aj could be rewritten as follows 

At = \Ai\(cos(6) + i • sin(<9)) = |^|exp(i • 6>), (58) 

where \A t \ is the modulus of A t and 6 = Arg(Ai) is the argument of A t . From Eq. (57), the following relation can be 
obtained 

^• = i-Log(|^|)-Arg(^). (59) 

The stability requires that Log(|/i;|) < 0. That means all eigenvalues of M 1 ™ 1 should lie within the unit circle the origin 
of which is located at the point (0, 0). The eigenvalue behavior of the matrix AT 1 with respect to the ^-parameter 
reflects the sensitivity of free parameters. 

3.2. Analysis of the sensitivities 

In this section, we consider the first- and second-order sensitivities of the eigenvalues with respect to the wave- 
number k and observe the real behaviors corresponding to v, n and c s in the L-MRT-LBM. Finally, in the complex 
plane, according to the distribution of the eigenvalues of the matrix A/™* and the stability condition, the influence of 
free parameters is investigated. 

3.2.1. The first- and second-order sensitivities with respect to the wave-number magnitude |k| 

Although there exist several papers that are focused on the analysis of dispersion and dissipation relations, as far 

as authors know, the investigations about dispersion and dissipation behaviors with respect to |k| still haven't been 

completely perfomed [3, 4, 7] for the MRT-LBM. 

From Eq. (15), the first-order derivatives_of acoustic modes (i.e. group velocity) 6c) ± (k) and shear modes 6i/(k) 

with respect to k = |k| are given as follows (k • u is the angle between k and u) 

Re[a>*(k)] = ±c s + |u|cos(k • u), 

7mK(k)] =-|k|(2gv + i7), 

Re[aj s K (k)] = |u|cos(k • u), 

7mK(k)] = -2|k|v. 

The second-order derivatives are given by 

{ M<(k)] = 0, 

/m[a&(k)] =-p^ v + J7 ), 

M<(k)] =0, 

7m[<(k)] =-2v. 
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Table 1: Illustrative parameters for the investigations of the original D2Q9 MRT-LBM (St = 1): e x denotes the unit vector (1,0) co-directional with 
x-axis. The symbol "\" indicates that the corresponding parameters are determined according to the cases. 



Index 




Se 


y 


n 


|u| 


k u 


u^ex 


k e x 





1.9 


1.64 


8.771929833e-3 


0.03658536587 


0.0 











1 


1.99 


1.64 


8.375209333e-4 


0.03658536587 


0.0 











2 


1.999 


1.64 


8.337503333e-5 


0.03658536587 


0.0 











3 


1.9999 


1.64 


8.333733333e-6 


0.03658536587 


0.0 











4 


1.99 


1.64 


8.375209333e-4 


0.03658536587 


0.1 


\ 


\ 


\ 


5 


1.9999 


1.64 


8.333733333e-6 


0.03658536587 


0.1 


\ 


\ 


\ 




0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 

K K 

Figure 1: Local information for derivatives of Im[aj]: (Left): Im[aj K ] ; (Righ): Im[aj KK ]. In the left subfigure, Green line denotes Im[aj s K ]; Red line 
denotes Im[aj^]. In the right subfigure, Green lines denote Im[aj^ K ]; Red line denotes Im[aj s KK ]. The parameters are from the index 1 of Tabel 1. 
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From Eqs. (60) and (61), we can detect the lattice sound speed c s (for |u| = or k • u = n/2) and the effective shear 
and bulk viscosities. Theoretically, the dispersion and dissipation relations are dependent on |k|, |u| and k • u for the 
fixed v and 77, and independent on the individual orientations of k and u. 

In order to observe the sensitivity of the hydrodynamic modes with respect to the wave-number /c, the locally- 
magnified profiles are given in Figs. 1 and 2. In the left sub figure of Fig. 1, the green profile stands for Im[oj s K ] and 
the red profile denotes Im[oj^]. It is noted that the acoustic dissipation has a significant change around k = 2. When k 
is larger than about 2.25, Im[oj^] grows very fastly. From the right sub-figure of Fig. 1, when k is smaller than about 
2.4, the shear viscosity v in the L-MRT-LBM will be smaller than the given v and becomes negative. This means 
that the real behavior of v and 77 depends on the wavenumber k. Meanwhile, the important fact is that the computed 
Im[oj^ K ] is varies from negative to positive values when k is larger that about 1.75. This phenomenon means that v + 77 
is changes from positive to negative values after k = 1.75 and the instability will appear for acoustic modes. In Fig. 2, 
the profiles of the left subfigure denote the numerical Re[oj K ]. The red lines are the exact Re[oj^]. It is clear that when 
the second term of Re[oj^] in Eq. (60), the red lines are the exact profiles of ±c s . From the subfigure, the computed 
sound speed c s is dependent on the wavenumber k under the condition of zero-mean flow. The green line indicates the 
exact Re[oj s K ]. It is clear that the computed Re[oj s K ] coincides with the exact Re[oj s K ] very well. From the right sub-figure 
in Fig. 2, we can see that computed Re[oj s K ] is exact for zero-mean flows. And the computed Re[oj^] is sensitive to 
the wave-number k. Especially, when the wave-number k is around 2, there exists some significant changes of the 
computed Re[oj^]. These results means the the numerical group velocity become very sensitive around k = 2. 




Figure 3: Local information for derivatives of Im\pj\\ (Left): Im[a> K ] ; (Righ): Im[aj KK ]. In the left subfigure, Green line denotes Im[aj s K ]; Red line 
denotes Im[aj^]. In the right subfigure, Green lines denote Im[a)^ K ]; Red line denotes Im[aj s KK ]. The parameters are from the index 2 of Tabel 1. 



In order to observe the impact of the influence of the small shear viscosity on the sensitivities, in Figs. 3~6, the 
profiles of the numerically-computed Im[a) K ], Im[a> KK ], Re[oj K ] and Re[oj KK ] are displayed for s v = 1.999 (indicated by 
the index 2 in Table 1) and s v = 1.9999 (indicated by the index 3 in Table 1). Clearly, the similar results are obtained 
compared with the results ( s v = 1.99, in the index 1 in Table 1). Now, we choose a smaller s v = 1.9 corresponding to 
the larger shear viscosity (indicated by the index in Table 1). From Figs. 7-8, it is observed that when k < 2, there 
does not exist any significant change for the sensitivities of the bulk viscosity v + 77 and the lattice sound speed c s . 
However, Fig. 7 shows that there exist some large deviations between the exact relations and the computed relations 
for Im[aj s K ] and Im[aj s KK ]. These deviations indicate that for the large shear viscosity v, the real v of the MRT-LBM 
deviates significantly from the exact v in a large range of wave numbers. 

In the above discussion, we focused on observing the sensitivity of oj with respect to k for the D2Q9 MRT-LBM. 
Now, we investigate the sensitivity of oj for the D3Q15 MRT-LBM. From Figs. 9~14, it is seen that there exist very 
significant changes for the first and second order derivatives of the acoustic and shear modes with respect to k. From 
Figs. 10, 12 and 14, it is observed that the computed c s has a significant dependence on the wave-number magnitudes, 
and the computed c s deviates significantly from the exact c s . From Figs. 9, 11 and 13, the computed Im[oj K ] and 
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Figure 4: Local information for derivatives of Re[aj]: (Left): Re[aj K ]; (Righ): Re[aj KK ]. In the left subfigure, Green lines denote Re[aj s K ]; Red line 
denotes Re[cjf]. In the right subfigure, Green lines denote Re[cj^ K ]; Red line denotes Re[co s KK \. The parameters are from the index 2 of Tabel 1. 



Table 2: Illustrative parameters for investigations of the D3Q19 MRT-LBM (St = 1): e x denotes the unit vector (1,0,0) codirectional with x-axis. 
The symbol "\ " indicates that the corresponding parameters are determined according to the cases. 



Index 




S e 


V 


1 


|u| 


k u 


u^ex 


k e x 





1.9 


1.6 


8.771929833e-3 


0.02111111118 


0.0 











1 


1.99 


1.6 


8.375209333e-4 


0.02777777778 


0.0 











2 


1.9999 


1.6 


8.333733333e-6 


0.02777777778 


0.0 











3 


1.99 


1.6 


8.375209333e-4 


0.02777777778 


0.1 


\ 


\ 


\ 


4 


1.9999 


1.6 


8.333733333e-6 


0.02777777778 


0.1 


\ 


\ 


\ 
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K K 

Figure 6: Local information for derivatives of Re[aj]: (Left): Re[aj K ]; (Righ): Re[aj KK ]. In the left subfigure, Green lines denote Re[aj s K ]; Red line 
denotes Re[cj^]. In the right subfigure, Green lines denote Re[a)^ K ]; Red line denotes Re[oj s KK ]. The parameters are from the index 3 of Tabel 1. 
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Figure 8: Local information for derivatives of Re[aj]: (Left): Re[aj K ]; (Righ): Re[aj KK ]. In the left subfigure, Green lines denote Re[aj s K ]; Red line 
denotes Re[cj^]. In the right subfigure, Green lines denote Re[a)^ K ]; Red line denotes Re[oj s KK ]. The parameters are from the index of Tabel 1. 




Figure 9: Local information for derivatives of Im\cj\\ (Left): Im[aj K ] ; (Righ): Im[aj KK ]. In the left subfigure, Green line denotes Im[aj s K ]; Red line 
denotes Im[aj^]. In the right subfigure, Green lines denote Im[aj* K ]; Red line denotes Im[oj s KK ]. The parameters are from the index 1 of Tabel 2. 



14 




Figure 10: Local information for derivatives of Re[aj]: (Left): Re[oj K ]; (Righ): Re[oj KK ]. In the left subfigure, Green lines denote Re[aj s K ]; Red line 
denotes Re[cj^]. In the right subfigure, Green lines denote Re[a)^ K ]; Red line denotes Re[oj s KK ]. The parameters are from the index 1 of Tabel 2. 




Figure 11: Local information for derivatives of Im[aj]: (Left): Imiaj^ ; (Righ): Im[cj KK ]. In the left subfigure, Green line denotes Im[aj s K ]; Red line 
denotes Im[aj^]. In the right subfigure, Green lines denote Im[oj^ K ]; Red line denotes Im[oj s KK ]. The parameters are from the index 2 of Tabel 2. 
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Figure 12: Local information for derivatives of Re[aj]: (Left): Re[oj K ]; (Righ): Re[oj KK ]. In the left subfigure, Green lines denote Re[aj s K ]; Red line 
denotes Re[cj^]. In the right subfigure, Green lines denote Re[a)^ K ]; Red line denotes Re[oj s KK ]. The parameters are from the index 2 of Tabel 2. 




Figure 13: Local information for derivatives of Im[aj]: (Left): Imiaj^ ; (Righ): Im[cj KK ]. In the left subfigure, Green line denotes Im[aj s K ]; Red line 
denotes Im[aj^]. In the right subfigure, Green lines denote Im[oj^ K ]; Red line denotes Im[oj s KK ]. The parameters are from the index of Tabel 2. 
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Figure 14: Local information for derivatives of Re[aj]: (Left): Re[oj K ]; (Righ): Re[oj KK ]. In the left subfigure, Green lines denote Re[aj s K ]; Red line 
denotes Re[aj^]. In the right subfigure, Green lines denote Re[a)^ K ]; Red line denotes Re[aj s KK ]. The parameters are from the index of Tabel 2. 



Im[co KK ] deviate from the exact Im[co K ] and Im[co KK ] . Numerically, Im[<jur] is very sensitive to the wave number k. 

From the current 2D and 3D results, using the given free relaxation parameters in the literatures, the numerical 
dispersion and dissipation relations are sensitive to the wave number k for u • k = 0. The behaviors of the numerical 
group velocity will become sensitive when the wave-number magnitude of |k| is close to 2. It is necessary to point out 
that if u • k ^ 0, we will obtain the most insensitive results. In order to save space, the investigations for |u| > and 
u • k ^ are not shown here. With respect to other free parameters, we will address these problems in complex plane 
in the next part. 

3.2.2. The behaviors of eigenvalues of M mn in the complex plane 

From Eq. (59), the linear stability of the MRT-LBM requires that the eigenvalues of M™ 1 should lie in a unit 
circle whose origin is located at the point (0,0). The investigations of sensitivities will address the factors that have 
a significant influence on the distribution of the eigenvalues in the complex planes. From Eq. (15), the distribution 
of the hydrodynamic modes are only dependent on the parameter set T = {k, u, k • u, v, rj}. That means the behavior 
of eigenvalues of the matrix M 1 ™ 1 should be dependent on the parameter set T. However, the researches dealing 
with the relations of dispersion and dissipation of the linearized MRT-LBM indicate that there exists some significant 
dependence on the orientations of k, u and free relaxation parameters in the matrix M™ 1 [7, 8]. The values of the free 
relaxation parameters are set according to the recommended values [7]. In order to implement the analysis, we give 
the following definations for k and u. For the 2D problems, k and u are defined by 

k x = |k|cos(^), k y = |k|sin(%); u = |u|cos(0 M ), v = |u|sin(0 M ). (62) 
The inner product of k and u is given by 

k • u = |k||u|cos(0* - 6 U ). 



k u = |k||u|cos(0* - o u ). (b3) 
For 3D problems, k and u are defined as follows 

k x = |k|cos(0jOsin(&), k y = |k|sin(0*)sin(&), K = |k|cos(&); (64) 
u = |u|cos(0 M )sin(^), v = |u|sin(0 M )sin(f M ), w = |u|cos(f M ). 3 

The inner product of k and u is given by 

k u = |k||u|(cos(0* - M )sin(&)sin(f M ) + cos(&)cos(f M )). (65) 

In Eq. (65), if let 6 k = 6 U , we have 

k-u = |k||u|cos(&-f M ). (66) 
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For the 2D problems, the discrete values of the variable |k| are set as a sequence of number {i • 8k}o<t<j where 
8k = n I J. First, the parameters given in Index 4 in Table 1 are used to implement our 2D investigations. Considering 
the wave-number vector k parallel to e x and k • e x is equal to or n, we consider k • u as a random parameter in 
the interval [0, n] . Ten uniform random values of u^e^ for each value of k are generated in [0, n] to investigate the 
influences of u^Sx on the distribution of the eigenvalues of the matrix M 1 ™ 1 . In order to address the influences of 
k • e x , we let u||e x ( u^e^ is equal to or n) , and ten uniform random values of k • e x for each value of k are generated 
in the range [0, tt]. In Fig. 15, the distribution of the eigenvalues of the matrix M™ 1 is given. From the figure, it is clear 
that k • e x for u||e x has more significant influence than u^e* for k||e x . When u is parallel to e x , The distribution of the 
eigenvalues of the matrix M 1 ™ 1 appears very random with respect to the random k • e x . In Fig. 16 , the distribution of 
the eigenvalues is figured for k • e x = n/4 and u^ex = tt/4. u • k is regarded as a random variable. In Fig. 17 , the 
distribution of the eigenvalues is displayed for k • e x = n/2 and vT^x = n/2. 

According to Figs. 15 - 17, we conclude that the distribution of the eigenvalues of the matrix M™^ is very sensitive 
to the angle between k and e x . Therefore, k and e x are expected to strongly govern the dissipation and dispersion 
relation of the linearized MRT-LBM. Furthermore, it is observed that the distribution of the eigenvalues appears 
symmetrical with respect to the real axis. This symmetry property is obvious because of the symmetry property of the 
wave-number vectors. The angle between the wave-number of the red color distribution and the wave-number of the 
blue color distribution is equal to n corresponding to each k. 




-2-1 1 2-2-1 1 2 

Re(l) Re(k) 

Figure 15: Distribution of the eigenvalues of the matrix M mrt : (Left) k||e x , and is regarded as a uniform random variable in the range [0, n]. 
For each value of k, ten random numbers for u^e* are generated in [0,n]. (Right) u||e x , and k • e x is regarded as a uniform random variable in 
the range [0, n]. For each value of k, ten random numbers for k • e x are generated in [0, n]. The parameters are from the index 4 of Tabel 1. (Blue 
Color) k = [k x , k y ]; (Red Color) k = [-k x , -k y ]. The circles in the subfigures stand for the unitary circles with the origin located at (0,0), which are 
also the stable regions of the linearized MRT-LBM. 

Here, we revisit the conclusion about the stability of the transverse modes and the longitude modes. It was 
indicated that the transverse modes is more stable than longitudinal modes and the sound waves propagating in the 
direction of the mean flow are quite unstable for the BGK-LBM [7]. This instability instability can be reduced in 
the original MRT-LBM [7]. In Fig. 18, when 0& and 6 U are regarded as the two random variables in the range [0, tt], 
the probability density function (PDF) profiles of of eigenvalue modules are given for the transverse modes and the 
longitudinal modes. From this figure, it is clear that for the original MRT-LBM, the PDF profiles of the eigenvalue 
modules of the matrix M™ 1 nearly has the same shape. And it is also observed that there exist some modules larger 
than 1 (refer to the vertical line x = 1 in the figure). That means the MRT-LBM is not always stable for any |k| e [0, tt] 
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-2-1 1 2-2-1 1 2 

Re(X) Re(k) 

Figure 16: Distribution of the eigenvalues of the matrix M 1 ™ 1 : (Left) k • e x = n/4; (Right) = n/4 . u • k is regarded as a uniform random 
variable in the range [0, n]. For each value of k, ten random numbers for u • k are generated in [0, n]. The parameters are from the index 4 of Tabel 
1. (Blue Color) k = [k x , k y ]; (Red Color) k = [-k x , -k y ]. The circles in the subfigures stand for the unitary circles with the origin located at (0,0), 
which are also the stable regions of the linearized MRT-LBM. 




-2-1 1 2-2-1 1 2 

Re(A.) Re(JL) 

Figure 17: Distribution of the eigenvalues of the matrix M™ 1 : (Left) k • e x = n/2; (Right) u"-~5 = tt/2 . u • k is regarded as a uniform random 
variable in the range [0, n]. For each value of k, ten random numbers for u • k are generated in [0, n]. The parameters are from the index 4 of Table 
1. (Blue Color) k = [k x , k y ]; (Red Color) k = [-k x , -k y ]. The circles in the sub figures stand for the unitary circles with the origin located at (0,0), 
which are also the stable regions of the linearized MRT-LBM. 
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\X\ 

Figure 18: Probability density functions of eigenvalue modules (In Eq. (62), % and 6 U are regarded as random variables in the rang [0,n] ) : (Red) 
k • e x = and u^e^ = 0; the difference of 6k and 6 U is equal to . (Blue) u • k = n/2 ; the difference of 6k and 6 U is equal to n/2. The parameters 
are chosen from Index 4 in Table 1 . 



and any small shear viscosity. These instability phenomena agree with the results of the second-order sensitivity 
analysis in Sec. 3.1.1. 

In Eq. (64), taking 6 k = U , we regard 0k and 6 U as random variables in the range [0, 2n]. Furthermore, we consider 
two cases: (a) assuming & = 0, £ u is a random variable in the range [0, 7r]; (b) assuming £ u = 0, & is a random variable 
in the range [0,7r]. In Fig. 19, the distribution of the eigenvalues of the matrix M 1 ™" 1 is plotted for any u • k. In Fig. 
19, the red points overlap with the blue points, meaning that for k = [k x , k y ] and k = [-k x , -k y ], similar eigenvalue 
behaviors exist. The random phenomenon in the right sub figure of Fig. 19 indicates that the eigenvalue distribution 
is very sensitive to the orientation of the wave number k. This phenomenon is identical with what was observed in 
the case of 2D linearized MRT-LBM. In Fig. 20, we consider two additional cases: (a) 0k = and & = 0, 6 U and g u are 
random in [0, 2n] and [0, n], repectively; (b) U = and £ u = 0, 0k and & are random in [0, 2n] and [0, n], repectively. 
From the figure, it is clear that we can obtain the same conclusion as that we have in Fig. 19 . 

From the investigations of eigenvalue distribution of the matrix M 10 ^ for the 2D and 3D linearized MRT-LBM, it 
is concluded that for the given u • k, the dispersion and dissipation behaviors are mainly determined by the orientation 
of k. Therefore, the acoustic modes and shear modes of the linearized MRT-LBM corresponding to the analytical 
hydrodynamic modes (15) are mainly sensitive to the wavenumber vector k. 

4. Simplified optimization strategies to determine the free relaxation parameters 

In this section, we will focus on simplifying the optimization strategy proposed in [8]. The simplification is based 
on the sensitivity analysis presented in Sec. 3 and the properties of the matrix B^ n \ 

An optimization strategy to obtain the optimal free relaxation parameters was proposed in [8]. This optimization 
strategy is well suited for the 2D/3D MRT-LBM. In order to handle low bulk viscosity flows, the researches in the 
present paper will be focused on dealing with the following minimization problem for 2D and 3D problems [8]: 

Min: £ 2d (H r )= f ||M^|||d^d(M)d^du, (67) 

Jo<e k <2n, 6t\k\<n, 0<e u <2n, \u\<u 

and 

Min : g M (Z r ) = f \\M ( :Y F d6 k dhdm)d6 u d£ u du, (68) 

JQ<6 k <27i, 0<&<tt, 6t\k\<n, 0<e u <2n, 0<£ u <n, \u\<u 
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Figure 19: Trajectories of the eigenvalues of the matrix M mt for the 3D linearized MRT-LBM. (Left) & = 0; (Right) = 0. The parameters are 
from the index 3 of Tabel 2. (Blue Color) k = [k x , k y ]; (Red Color) k = [-k x , -k y ]. The circles in the subfigures stand for the unitary circles with 
the origin located at (0,0), which are also the stable regions of the linearized MRT-LBM. 
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Figure 20: Trajectories of the eigenvalues of the matrix M™" 1 for the 3D linearized MRT-LBM. (Left) 6k = and & = 0, 9 U and £ u are random in 
[0, 2n] and [0, n], repectively; (Right) 6 U = and £ u = 0, % and & are random in [0, 27r] and [0, 7r], repectively. The parameters are from the index 
3 of Tabel 2. (Blue Color) k = [k x , k y ]; (Red Color) k = [-k x , -k y ]. The circles in the subfigures stand for the unitary circles with the origin located 
at (0,0), which are also the stable regions of the linearized MRT-LBM. 
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where u indicates |u|, H r denotes the free relaxation parameter set ( H r c H), || • \\p indicates the matrix Frobenius norm 
and the matrix M e is dependent on |u|, |k|, St and free relaxation parameters, which is defined by 



M (n) = g(n) _£ + St 2^ 



In Eq. (69), & n) and & are defined by 



The matrix & in Eq. (69) is given as follows for the 2D/3D L-MRT-LBM 











-\k x 2 cr e 

^k X kyCJ~ g 





3 e 



(69) 
(70) 

(71) 



and 



-\k y k x {cr v + 2cr e ) 



-\k y k x (cr v + 2cr e ) 
-\k y 2 {cr v + 2cr e ) 



k z k x (cr v + 2cr e ) - gk z k y (cr v + 2cr e ) 



-\k z k x (cr v + 2cr e ) 
-\k z ky{cr v + 2cr e ) 



- l -k 



'(O-y + 2(T Z ) 



(72) 



which are corresponding to the exact total bulk viscosity dissipation coefficient matrix of St (refer to Appendix B ). 

According to the analysis displayed in Sec. 3.2, the eigenvalues are mainly sensitive to the wave number k. 
Therefore, taking u to be parallel to e x and considering the most sensitive factors, we can simplify Eq.(67) and Eq. 
(68) as follows 



Min 



*>-/ 

Jo<_ 



||Al< n) ||ld0*d(M)du, 



0<6 k <2n, St\k\<7:, |u-e x |<M 



and 



Min : 



J0<> 



\\M[ n) \\ 2 F de k d{ k d(6tk)du, 



(73) 



(74) 



J0<e k <2n, 0<£ k <n, 6t\k\<n, |u-e x |<w 

where u = |u • e x |. More specifically, the above simplified formula Eq. (74) is very effective for obtaining the free 
relaxation parameters of the 3D problems, which reduce the integral variables. 

In Table 3, illustrative examples are given for both 2D and 3D problems. The truncated error of Eq. (38) is set 
to equal 4 for handling the acoustic problems with the low bulk viscosity and for guaranteeing stability [8]. In order 
to investigate the dispersion and dissipation relations for the optimal MRT-LBM, the profiles of the transverse modes 
and the longitudinal modes are shown in Figs. 21-24. From these figures, it is clear that the simplified optimization 
strategies (73)~ (74) can make the MRT-LBM schemes endowed with the low-dispersion low-dissipation properties. 
From the viewpoint of stability, the optimal MRT-LBM appears more stable than the original MRT-LBM [7, 8]. 
Especially, for the optimal MRT-LBM, the lower bulk dissipation parameter can be chosen for aeroacoustic problems. 
In the original MRT-LBM, because of sustaining the stability of the MRT-LBM schemes, the low bulk viscosity is 
forbidden or impossible [7, 8]. From Figs. 21 - 23, it is observed that for the optimal MRT-LBM, the acoustic modes 
with u perpendicular to k are more stable than the acoustic modes with u parallel to k, and the shear modes with u 
perpendicular to k are more accurate than the shear modes with u parallel to k. 



5. Numerical examples 

In this section, we perform some numerical simulations to validate the optimal MRT-LBM on 2D/3D benchmark 
problems . All simulations are performed using an ad hoc version of PalaBos [20]. For the current optimization 
procedures, the truncated error term is only up to the fourth order so that we can improve the dissipation relation 
errors of the MRT-LBM and the isotropic errors [8]. In this part, we only focus on the case with the fourth-order 
truncated error terms. If you want to further improve both the dispersion and dissipation errors, the higher-order 
truncated error terms are needed [8]. 
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Table 3: Optimized free relaxation parameters corresponding to the specified s v and s e for the MRT-LBM under the condition u||e x (u^e^ = 0) . 
The upper bound uq of |u| is equal to 0.2. The symbol "\ " indicates that the corresponding fields are inapplicable. For D2Q9, e x = (1, 0). For 



D3Q19, e x = (1, 0, 0). The truncated error in Eq. (38) is equal to St 4 . 

Model Index s v s e s e s q s n s t 

D2Q9 1.99960008 1.99960008 1.997623852 1.999448768 \ \ 

D3Q19 1 1.99960008 1.99960008 1.994155318 1.998982400 1.992025433 1.999363601 

D2Q9 2 1.886792453 1.886792453 1.470217043 1.842711354 \ \ 

D3Q19 3 1.999960001 1.9999800002 1.9993910550 1.999916523 1.9996934332 1.999886136 

D2Q9 4 1.995389843 1.9953898430 1.9729857787 1.9875516898 \ \ 

5 1.999996000 1.9999960008 1.9999762501 1.999994487 \ \ 

6 1.9999988 1.9999988 1.9999928163 1.9999992651 \ \ 
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(a) u||k (b) u_Lk 

Figure 22: The dispersion relations for the optimized 2D MRT-LBM: (a) = 0; (b) = 0, k • e x = n/2. The relaxation parameters are 
given by Index in Table 3. The magnitude of u is equal to 0.2. (Green line) Exact shear modes Re[c*/]; (Red line) Exact acoustic modes Refw*]. 
co ±,e and aj s,e denote the exact acoustic modes and shear modes respectively. and cj s denote the numerical acoustic modes and shear modes 
respectively. 




24 



(a) u||k 



(b) u_Lk 



Figure 24: The dispersion relations for the optimized 3D MRT-LBM: (a) u • e x = 0; (b) u • e x = 0, 6k = U = and & = n/2 in Eq. (64). The 
relaxation parameters are given by Index 1 in Table 3. The magnitude of u is equal to 0.2. (Green line) Exact shear modes Re[6/]; (Red line) Exact 
acoustic modes Refw*]. a) ±,e and cj s,e denote the exact acoustic modes and shear modes respectively. a) ± and cj s denote the numerical acoustic 
modes and shear modes respectively. 



5.7. A 2D acoustic point source 

In this part, we consider propagation of the 2D acoustic wave generated by a point source. These investigations 
only focus on comparing the optimal MRT-LBM with the BGK-LBM and the quasi-incompressible BGK-LBM. The 
time-dependent acoustic point source is set as 

p(-,0 = 1 +p/sin(27r/r • t) 

u(-, t) =0 (75) 
v(-,0 =0 

where p/ = 0.01. The compuational domain is £1 = [0, L] 2 (L = 1). The number of the lattice points is equal to 
200 x 200. The period T is equal to 4 • Sx (Sx - 1/200). The relaxation parameters in the standard LBM are set by 
Index 5 in Table 3. In Fig. 25, the contour lines for three different methods are shown based on the same contour 
level values. It is very clear that the contour line result obtained by the optimal MRT-LBM is the best one. The 
BGK-LBM has nearly the same results as the quasi-incompressible BGK-LBM. This results are similar to that in 
former researches [8]. Here, we point out that, with respect to the optimal MRT-LBM, the contour line result of the 
acoustic point source in our former research is better than the result of the optimal MRT-LBM in the current paper. 
It is needed to be clarified that in the current optimization strategy, the mean flow influence is considered, but in 
the former research, the optimization strategy for the kind of problem was implemented with the zero mean flows. 
The optimization strategies in this paper appear in a more general form. From the results of this problem, the well 
performances of the optimization strategies are further demonstrated. 

5.2. The 2D acoustic pulse source 

We now consider the 2D acoustic pulse source for assessing the simplified 2D optimization strategy. Assuming 
that the viscosity effect is negligible on acoustic waves, the acoustic pulse problem possesses an analytical solution 
[8, 21]. The initial profile is given as follows 

po =l+p/ 

u = U (76) 
v =0 
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Figure 25: Density contour lines for the BGK-LBM, the quasi-incompressible BGK-LBM and the optimal MRT-LBM at t = 0.7. The viscosity 
v= 1.66667 x 10" 5 . 
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where p/, e, a, r and Uq are defined by 

pr = 6exp(-a • r\ e = 10" 3 , a = ln(2)/Z? 2 , r = yj(x - x ) 2 + (y - y ) 2 , U = O.Ol. (77) 
The parameter b in Eq. (77) is equal to 0.03. The exact solution of pr (if (x§,y§) - (0, 0)) is given by [21] 

P'(x,y,z) = ^j Q exp|-£jcos(c s ^J (fi7)fdf (78) 



where r] = ^(x- U^t) 2 +y 2 + z 2 and Jq(-) is the Bessel zero-order function of the first kind. The computational 
domain Q = [0, L] 2 (L = 1). The relaxation parameters are determined by Index 5 in Table 3. In order to validate 
the precision of the optimal MRT-LBM, the mesh convergence and the time convergence are calculated. The adopted 
lattice resolutions are 100 2 , 200 2 , 300 2 and 400 2 . The evolution of the acoustic pulse is restricted in the computational 
domain in order to avoid the effects of boundaries. From here on, the L 2 relative error is defined by 



E L 2(t) 



Z£T(P'th(Xf, - P'numfo, t)) 2 

where p'thfe and p/ n um(x/, t) denote the theoretical and numerical solutions at the lattice nodes x,-, respectively. 
First, in Fig. 26, the density profiles for different LBS are given at t = 0.4. It is clear that the results obtained by 
the original MRT-LBM are less satisfactory because of the large bulk viscosity. In order to show the quantitative 
comparisons, the L 2 relative errors are given in Figs. 27-28 with respect to the mesh resolution and the evolution 
time, respectively. From Fig. 27, it is clear that the results obtained by the optimal MRT-LBM have the lowest 
L 2 relative errors and the best convergence property. Further, some super-convergence properties are exhibited for 
the optimal MRT-LBM. Let us observe Fig. 28 and we can conclude that when the mesh resolutions are fixed, the 
original MRT-LBM has the best precision. With respect to the evolution time, three methods have nearly the same 
time convergence order. 

5.3. The 3D acoustic spherical pulse source 

In this part, we continue to investigate the 3D acoustic spherical pulse source in order to obtain a further validation 
for the simplified 3D optimization strategies. Assuming the viscosity influence is ignored, the acoustic pulse problem 
possesses the analytical solution [22]. The initial profile is given as follows 

(80) 

wo =0 

where p/, e, a, r and Uq are defined by 

p/ = eexp(-a • r 2 ), e = 10" 3 , a = ln(2)/Z? 2 , r = ^(jc - x ) 2 + (y - y ) 2 + (z - z ) 2 , U = 0.01. (81) 
The parameter b in Eq. (81) is equal to 0.03. The exact solution of pr (if (xq, yo, zo) = (0, 0, 0)) is given by [22] 

P'(x,y 9 z)=T-^= f exp(-|^)cos(c^)^^^ (82) 

where rj = ^/(x - Uot) 2 + y 2 + z 2 . The computational domain Q = [0, L] 3 (L = 1). The simulation results of two 
cases are shown. The relaxation parameters are determined by Index 3 in Table 3. Because of the very low viscosity, 
the corresponding viscosity dissipation can be ignored. The adopted lattice resolutions are 100 2 , 200 2 and 300 2 . In 
Fig. 29, the density profiles are shown for different LBS. It is observed that the profiles obtained by the original 
MRT-LBM have the worst precision. Quantitatively, the L 2 relative errors with respect to the mesh resolutions and 
the evolution time are presented in Figs. 30-31. The conclusions are similar to that of the 2D problems. The optimal 
MRT-LBM has the best precision with respect to the mesh resolutions and the evolution time. At the same time, some 
super-convergence property of the optimal MRT-LBM is observed. 

In all, from the numerical investigations, the optimal MRT-LBM not only improves the computational precision 
for acoustic problems, but also present some super-convergence properties compared with the BGK-LBM. 



Po 


= 1 + p/ 
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= U 
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Figure 26: Fluctuation density profiles of 1000p/ for three different 2D LBS along the line y=0.5 at t=0.4: (a) Lattice number 100 2 ; (b) Lattice 
number 200 2 ; (c) Lattice number 300 2 ; (d) Lattice number 400 2 . 
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Figure 27: The L relative errors for three different 2D LBS with respect to the mesh scale 6x: BGK-LBM, Optimal MRT-LBM, Original MRT- 
LBM. 
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10" 1 




(c) Lattice 300 2 (d)Lattice 200 2 



Figure 28: The L 2 relative errors for three different 2D LBS with respect to the evolution time t: BGK-LBM, Optimal MRT-LBM, Original 
MRT-LBM. Lattice number: (a) 100 2 ; (b) 200 2 ; (c) 300 2 ; (d) 400 2 . 
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Figure 29: Fluctuation density profiles of 1000p/ for three different 3D LBS along the line y=z=0.5 at t=0.4: (a) Lattice number 100 3 ; (b) Lattice 
number 200 3 ; (c) Lattice number 300 3 . 
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Mesh scale 8 x 



(c) at t=0.5 

Figure 30: The L? relative errors for three different 3D LBS with respect to the mesh scale 6x: BGK-LBM, Optimal MRT-LBM, Original MRT- 
LBM. 
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(a) Lattice 100 3 (b)Lattice 200 3 
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(c) Lattice 300 3 



Figure 31: The L 2 relative errors for three different 3D LBS with respect to the evolution time t: BGK-LBM, Optimal MRT-LBM, Original 
MRT-LBM. Lattice number: (a) 100 3 ; (b) 200 3 ; (c) 300 3 . 
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5.4. The sound generation by a turbulent airflow over a moderately deep rectangular cavity 

The problem is to study aerodynamic sound generation by a turbulent air flow over a moderately deep cavity. 
Despite the real flow is expected to be 3D and turbulent, it is known that dominant acoustic tonal modes are generated 
by nearly 2D structures which can be captured thanks to 2D simulations. Therefore, 2D simulations may be useful in 
predicting main frequencies of the radiated acoustic field. This benchmark problem is taken from the "Fourth Com- 
putational Aeroacoustics (CAA) Workshop on Benchmark Problems". The ratio of the cavity width and the depth is 
0.3124. In Fig. 32, the schematic view of the computational domain is given along with the corresponding physical 
scales. The inflow at the far left boundary was set as uniform with a velocity magnitude of 50 m/s along the horizon- 
tal direction. The outflow boundary condition and free- slip boundary condition are applied at the outlet and the top 
boundary, respectively. In order to avoid the influence of the sound wave reflection from the left boundary, the outlet 
boundary and the top boundary, buffer (sponge) regions are used. The buffer region thickness at the inlet and outlet 
is equal to one-twentieth of the domain length along x-direction. The buffer region thickness for the top boundary is 
one-twentieth of the rectangle region. The lattice resolution along the thickness of the cavity tongue B is equal to 50. 
The optimal relaxation parameters are set by Index 4 in Table 3. It is necessary to point out that for the current con- 
figuration, the computation by the BGK-LBM diverged after about 1.4 x 10 5 time steps. For the current investigation, 
the computation by the optimal MRT-LBM was implemented for 6.48 x 10 5 time steps and the optimal MRT-LBM 
was still stable. Figs. 33-35 show the computed Fourier spectra of sound pressures with the measured sound pressure 
spectra and the CFD (TIDAL, Time Iterative Density/pressure based Algorithm) sound pressure spectra. The tonal 
frequencies and amplitudes of the computed spectra coincide with those of the experimental spectra very well. It is 
clear the obtained spectra by the optimal MRT-LBM is better than that by TIDAL. 
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Figure 32: Schematic computational geometry: A = 25.42mm, D = 7.94mm, B = 3.18mm. 



5.5. Application to the acoustic scattering from a cylinder solid body 

The acoustic scattering from a cylinder solid body was presented in the second CAA workshop on benchmark 
problems [23]. The objective of this problem focused on addressing the propeller generated sound field scattered 
off the fuselage of an aircraft. By this benchmark, we address the performance comparisons among the optimal 
MRT-LBM, the DRP schemes and the high-order compact schemes for solving the linear acoustic problems (the 
DRP schemes and the high-order compact schemes are used to solve the linearized Euler equations). The schematic 
geometry is given in Fig. 36. The left big black disc indicates the cylinder and the right small black disc denotes the 
acoustic source. The two-dimensional cylinder radius R is equal to D/2 (D = 1 for current researches). The centre of 
the cylinder is located at the origin. At the initial time (t - 0), the initial profile of the acoustic source is defined by 



p(x,y,0) = poo + £exp 



-2 



(x-x s Df+y 2 



(83) 
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Figure 33: Comparison of sound pressure spectra at the centers of the cavity left wall 
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Figure 34: Comparison of sound pressure spectra at the centers of the cavity floor 
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Figure 35: Comparison of sound pressure spectra at the centers of the cavity right wall 



where w = 0.2, 6 = c 2 x 10" 4 , x s = 4 and = c 2 poo- The low 6 is to ensure a linear response and guarantee the 
linear analysis theory. The problem was propsed as the second problem in Category 1 in [23], but the given analytic 
solution is incorrect [24]. The correct analytical solution was given as follows [24, 25] 

p e (x, y, t) = poo + Re | ^ (Ai(x, y, co) + A r (x, y, co)) ^exp(-k^0d^ j . (84) 

The incident wave amplitude is 

At(x, y, co) = ^-exp{-co 2 /(4b)}J (cor s ) (85) 
lb 

where r s = ^(x- x s D) 2 + y 2 , b = ln2/w 2 and Jo is the Bessel function of order zero. And the corresponding reflected 
wave portion is given by 

oo 

A r (x, y, co) = A r (r, co) = ^ C k (co)H^\rco)cos(kco), (86) 



k=0 



where is the kth order Hankel function of the first kind and 



C k (co) = ^exp{-co 2 /(4b)} f J^corj '" ~ s ' cos(fcfl)dfl, (87) 

Jo 



£k f °° T , x s Dcos(6) 



where e = l,e k = 2 for k ± 0, r = 1/2 and r So = ^0.25 + x 2 D 2 - x s Dcos(6). 

In order to compare numerical results with the exact solution, the computational time is non-dimensionalized as 
tc s /D . The pressure samples are recorded at three different points A(r = 5D, 6 = 7r/2), A(r = 5D, 6 - (3/4)n) and 
C(r = 5D, 6 = 7r). The pressure signals are sampled between 5 < tc s /D < 10. The acoustic pressure p osc is non- 
dimensionalized as p' - p osc / [(pcoC 2 )S^ . For the current computations, the relaxation parameters are determined by 
Index 6 in Table 3 and the lattice resolution along the diameter of the cylinder is equal to 50. The corresponding values 
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of the shear and bulk viscosity are very small, so that the dissipation influence on the acoustic waves can be ignored. 
In Fig. 37, the comparisons between the optimal MRT-LBM, the DRP schemes and the high-order compact schemes 
are presented. For the schemes of the DRP schemes and the high-order compact schemes, the temporal integration 
is performed using the 4-6 low dispersion and dissipation Runge-Kutta (LDDRK) optimized scheme [26, 27]. From 
the results, it is clear that the best results are obtained by the optimal MRT-LBM. This conclusion guarantees that 
the optimal MRT-LBM is better than the fourth-order DRP scheme and sixth/fourth-order compact schemes of the 
linearized Euler equations. In Fig. 38, the comparisons between the optimal MRT-LBM and the exact solution are 
shown. From the results, it is observed that the computational error at the points B and C is larger than that at the point 
A. The observed deviations from the exact solution are mainly from the non-linear influence of the lattice Boltzmann 
schemes, because the exact solution is obtained by the linearized Euler equations. The another possible reason is 
attributed to the step approximation of the the cylinder boundary. However, the results of the optimal MRT-LBM are 
the best ones. In Fig. 39, the acoustic wave patterns are shown at tc s /D = 5.77 and tc s /D = 8.66. Three wave fronts 
are observed. The farthest wave front from the cylinder is generated by the initial condition. The next front is the 
wave reflected off the right surface of the cylinder. The wave front closest to the cylinder is produced by the collision 
and merging of the two parts of the initial wave front on the left side of the cylinder. 




Figure 36: Schematic figure of the acoustic scattering problem. 



6. Conclusion 



In this paper, the researches are focused on the sensitivity analysis and the determination of the relaxation param- 
eters in the MRT-LBM. By the investigations of the sensitivity, it is concluded that the numerical behaviors of the 
dispersion and dissipation are mainly dependent on the variation of the wave number k. From the behaviors of the 
eigenvalues of M mrt in the complex planes, it is observed that the distribution of the eigenvalues has the certain shapes 
if the wave number k is fixed and the velocity vector is regarded as a random variable. According to the sensitivity 
analysis, the integral spaces of the optimization strategies are simplified for 2D/3D problems. Then, by von Neu- 
mann analysis, the numerical performances of the optimization strategies are validated. Finally, by the benchmark 
acoustic problems and considering the ignorable small viscosity in the LBS, the optimal MRT-LBM is compared 
with the BGK-LBM and the used high-order compact/DRP schemes . The numerical results demonstrate that the 
optimal MRT-LBM can be used to simulate the acoustic problems very well. Meanwhile, the optimal MRT-LBM 
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Figure 37: Comparison between computed (the optimal MRT-LBM, the DRP schemes and the high-order compact schemes) and exact solutions 
of two-dimensional scattering problem at the point A . (a) Compared with the fourth-order DRP scheme (seven-point stencil), (b) Compared with 
the sixth-order compact scheme (three-point stencil), (c) Compared with the fourth-order compact scheme (three-point stencil). Line ( — ): Exact 
solution. Dashed line ( ): optimal MRT-LBM solution. Dash-dot line (- • -): solutions of the high-order/DRP schemes. 
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Figure 38: Comparison between computed (optimal MRT-LBM) and exact solutions of two-dimensional scattering problem at the three points . 
(a) The time series of p' at the point A. (b) The time series of p' at the point B. (c) The time series of p' at the point C. Line ( — ): Exact solution. 
Dashed line ( ): optimal MRT-LBM solutions. 




Figure 39: The acoustic wave patterns by the optimal MRT-LBM. (Left) tc s /D = 5.77; (Right) tc s /D = 8.66. 
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performs better than the BGK-LBM and the used high-order compact/DRP schemes numerically. Especially, some 
super-convergence properties are observed. 
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Appendix A. Transformation matrices and matrices *P of MRT-LBM 

Appendix A.l. For two dimensional MRT-LBM with 9 discrete velocities 

The transformation matrice for MRT-LBM with 9 discrete velocities is described by 
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(A.1) 



For the weakly-compressible MRT-LBM, the matrix *F is given by 
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For the MRT-LBM similar to an incompressible LBS, the matrix *F is given by 
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Appendix A.2. For three dimensional MRT-LBM with 15 discrete velocities 

The transformation matrice for MRT-LBM with 15 discrete velocities is described by 
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Appendix B. Coefficient Matrices 

In this section, the dissipation coefficient matrix will be divided into two parts: exact dissipation coefficient matrix 
(St • D)\ the part (St • D) of non-Galilean invariance related with mean flows. 

Appendix B.l. Coefficient matrices ofL-NSE with the first-order of St for 2D MRT-LBM 
The exact part D of the dissipation coefficient matrix is given by 





D 



3 kxkyCT e 



3 kxkyO~e 3 ^x @~v 3 ky @~v 3 ky @~e 
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(B.l) 



The part D of non-Galilean invariance is expressed by 

D[l, 1] = D[l, 2] = D[l, 3] = (B.2) 
5[2, 1] = \k x 2 Ucr v + \k y 2 Ucr v + \k x k y Vcr e + \k x 2 Ucr e 

D[2, 2] = k x 2 V 2 cr e + 2k x 2 U 2 cr e + k x Vk y Ucr e - k x 2 V 2 cr v + 2k x 2 U 2 cr v + ^VJkyt/ov + fc/vV y 
D[2, 3] = 2k x k y V 2 cr e + k x 2 UVcr e + k x U 2 k y cr e + 2k x U 2 k y cr v - 2k x k y V 2 cr v + 3ky 2 VU<r v - k x 2 UVcr v 

(B.3) 

D[3, 1] = ^/Vcr, + \k y k x Ucr e + |£/Vcr y + ^ 2 Vo- y 

5[3, 2] = k y 2 VUcr e + k x k y V 2 cr e + 2k x U 2 k y cr e + 3£ X 2 £/Vcr y + 2k x k y V 2 (T v - 2k x U 2 k y cr v - k y 2 VUcr v 
D[3, 3] = & v 2 t/V e + k x Vk y Ucr e + 2ky 2 V 2 (r e - k y 2 U 2 cr v + 4^V^t/o- y + 2£ v 2 VV y + k x 2 U 2 cr v 

Appendix B.2. For three dimensional MRT-LBM with 19 discrete velocities 



D = 



-^ y |k| 2 -^ x Vv + 2cr,) 



~\k y k x (cr y + 2cr,) -|o- y |k| 2 - ^/(cr y + 2cr,) ~\k z k y (cr v + 2cr,) 

~\k z k x (cr v + 2o- e ) -^ z fcy (o~ v + 2cr e ) -|cr y |k| 2 - \k 2 (cr y + 2cr z ) 



(B.4) 



where the wave-number vector k = (k x ,k y , k z ) and |k| = ^k 2 x + & 2 + & 2 . The viscosity is defined by v = |<x y and the 
bulk viscosity is defined by r] = \o~ e . Then, the effective bulk viscosity is given by n + \v. 



The part D of non-Galilean invariance is given by 

5[1, 1] = 2] = D[l, 3] = 4] = 



(B.5) 
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D[2, 1] = \k y 2 Ucr v + \k 2 Ucr y + \k x Vk y cr y + \k x Wk z cr y + \k 2 Ucr y + \k x Wk z <r e + \k x Vky<r e + \k 2 Ucr e 
D[2, 2] = fk x Wk z Ucr v + 2ik z Vik y Wo- v + jfc z ¥V v + fk x Vk y Ucr y + fik/t/Vy - fik/vVy + Jfc^VVy- 

ffc/wVy + fkW^tftre + f fc/t/V, + \k 2 W 2 a e + f^V^t/^ + ffc/yV, 
D[2, 3] = \k y k x U 2 o- y + 3*y Wfc z t/cr v - f k x k y V 2 (r v - lk x W 2 k y cr y - \k 2 UV(T y + 3^ 2 yt/o- y - \k x Wk z V(T y + 

\k x U 2 k y (T e + IkxWPkyO-e + IkxWkzVo-e + §£ x 2 t/ycr, + fyt^vV, 
5[2, 4] = -§ k x 2 UWcr v + 3£ z 2 Wo- y - 2/3^V^Wo" y - \k x k z W 2 a v + 3k y Uk z Vcr v + \k x U 2 k z (T y - \k x V 2 k z (T y + 

lk 2 UWo- e + \k x U 2 k z o- e + fy^V^Wo-, + f Jk^wV* + f^V 2 ^ 
£>[3, 1] = l^/Vcr,, + \kyWk z (T e + \k y Uk x <j e + |& z 2 y<x y + \kyWk z cr y + \k y Uk x cr y + |fcy 2 Vcr v + ^ x 2 Vcr v 
D[3, 2] = p x U 2 k y cr e + |^W 2 ^o- e + ffc/yt/o^ + f^W^t/cr* + fJkjJfcy VV* - |^W 2 ^cr y - \k y k x U 2 (T y + 

3k x 2 UVcr v + 3^W^Vo- y - lk y 2 VU(T v + ^^VV y - |^W^t/o- y 
D[3, 3] = f yk/wV, + \k y 2 V 2 (T e + + \kyWk z Vo- e + \k y 2 U 2 (T e + fk x Vk y Ucr y + 2Jfc JC WJfc z t/o- v + 

k x 2 U 2 cr v + ^¥V V -U*t- lk y 2 W 2 (T v - \ky 2 U 2 cr y + |^ 2 y 2 cr y + ^ik^Wo-v 
D[3, 4] = 3k z Uk x Vcr v - \k y U 2 k z (Ty - lk y 2 VWcr y - \k y k z W 2 (r y - \k y Uk x W(T y + \kyV 2 k z (T y + 3£ z 2 Wcr v + 

\k y U 2 k z cr e + f /^WV, + \ky\ 2 k z cr e + f^t/^Wo^ + |^ 2 VWo-, 
D[4, 1] = |^ z 2 Wcr e + ^k z Vk y cr e + \k z Uk x <j e + |& z 2 W<x y + \ky 2 W<r y + \k z Uk x cr y + |^ 2 \y<r y + ^k z Vk y cr y 
D[4, 2] = lk x k z W 2 o- e + f^y^t/cr, + \k x U 2 k z o- e + f^y 2 ^ + §£ z 2 Wcr, + 3& x 2 t/Wo- y - §£ z 2 Wcr v + 

3^y^Wo- y + \k x k z W 2 cr y - \k x \ 2 k z cr y - \k x U 2 k z cr y - lk y Uk z Vcr y 
D[4, 3] = lk z Uk x Vcr e + §£ z 2 Wcr, + f k y V 2 k z cr e + fy^WV, + \k y U 2 k z (T e - lk 2 WVcr y + 3^ 2 yWcr y - 

\k y U 2 k z (T y + 3kyUk x Wo- Y - \kyV 2 k z o- y + Ikyk^o-y - \k z Uk x Vo- y 
D[4, 4] = f £ z 2 W 2 cr v + fk x Wk z Uo- v + £ v 2 yV y + 2k x Vk y Ucr v - \k 2 U 2 cr y + k 2 U 2 cr y + f Jk^Jky Wcr v - 

f ^ z 2 y 2 o- y + |^^ z t/o- e + |^ z 2 t/ 2 o-, + f £ Z 2 WV, + |/: z 2 y 2 o- e + |jkyWik z y^ 

(B.6) 
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